POC MICROSOFT-Ibema

Sumário Executivo

A Ibema, juntamente com a Microsoft, solicitaram a BlueShift que avaliasse a possibilidade de melhorar as estratégias de go-to-market da Ibema por meio da aplicação de técnicas de inteligência artificial e de machine learning aos dados de volumes de vendas e de preços da Ibema.

A fim de atender a essa demanda, a BlueShift criou modelos de forecasting de volumes e de preços líquidos para a Ibema por meio de técnicas tanto consagradas quanto de ponta, a fim de apoiar o processo decisório da empresa no que diz respeito à otimização de esforços e de alocação de recursos para regiões e sites de maior interesse, ou seja, aqueles com melhor prognóstico de faturamento.

A aplicação das técnicas mencionadas anteriormente para os dados do estado de São Paulo, bobinas do site TVO, permite concluir que os modelos preditivos desenvolvidos apresentaram muito bom desempenho e que, portanto, podem ser usados para prever acuradamente tanto o preço líquido quanto o volume de vendas em um horizonte de até 14 dias. Essas previsões, se estendidas para outras regiões, sites e tipos de embalagem, certamente poderão apoiar o processo decisório da Ibema no que diz respeito ao planejamento e controle da produção, bem como em relação às estratégias de marketing.

Introdução

O objetivo desta PoC é criar modelos de forecasting de volumes e de preços líquidos para a Ibema por meio de técnicas tanto consagradas quanto de ponta, a fim de apoiar o processo decisório da empresa no que diz respeito à otimização de esforços e de recursos para regiões e sites de maior interesse.

Será usada como métrica de interesse de uma região/site a previsão de faturamento líquido - preço líquido x volume - nos próximos 14 dias.

Ingestão e Preparação dos Dados

Os dados recebidos da Ibema em formato xlsx foram ingeridos e, depois de alguns esforços de análise exploratória, receberam o seguinte tratamento:

• Variáveis tidas por categóricas: cod_site, cod_produto, cod_uf, cod_tipo_embalagem, cod_gramatura, cod_configuracao

• Variáveis tomadas como estritamente positivas: num_comprimento, num_largura, qtd_volume, vlr_faturamento_bruto, vlr_preco_liquido

Variáveis não utilizadas: vlr_despesa_financeira, vlr_faturamento_bruto, vlr_comissao, vlr_despesa_venda, vlr_frete_unitario, nom_cidade, vlr_faturamento_sem_imposto, vlr_impostos, vlr_taxa_nota, vlr_faturamento_sem_ipi, vlr_frete

A seguir, mostram-se algumas estatísticas descritivas do dataset resultante dessa etapa, considerando:

• UF: ‘SP’
• SITE: ‘TVO’
• EMBALAGEM: ‘BOBINA’

UF        <- 'SP'
SITE      <- 'TVO'
EMBALAGEM <- 'BOBINA'

ibema <- read_delim(
  "00_data/Dados_BlueShift_Total.csv", ";",
  escape_double = FALSE,
  col_types = cols(
    VLR_PRECO_LIQUIDO           = col_number(),
    VLR_FATURAMENTO_BRUTO       = col_number(),
    VLR_FATURAMENTO_SEM_IPI     = col_number(),
    VLR_FATURAMENTO_SEM_IMPOSTO = col_number(),
    VLR_TAXA_NOTA               = col_number(),
    QTD_VOLUME                  = col_number()
  ),
  trim_ws = TRUE
) %>% 
  janitor::clean_names() %>% 
  mutate(
    cod_gramatura      = as.factor(cod_gramatura),
    cod_site           = as.factor(cod_site),
    cod_uf             = as.factor(cod_uf),
    cod_tipo_embalagem = as.factor(cod_tipo_embalagem),
    cod_configuracao   = as.factor(cod_configuracao),
    cod_produto        = as.factor(cod_produto),
    dat_movimento      = as.Date(dat_movimento),
    cod_produto        = fct_lump_prop(cod_produto, prop = 0.01),
    cod_configuracao   = fct_lump_prop(cod_configuracao, prop = 0.01),
    cod_uf             = fct_lump_prop(cod_uf, prop = 0.01),
    cod_gramatura      = fct_lump_prop(cod_gramatura, prop = 0.01)
  ) %>% 
  filter(
    num_comprimento > 0, num_largura > 0,  qtd_volume > 0, 
    vlr_faturamento_bruto > 0, vlr_preco_liquido > 0
  ) %>% 
  select(
    -vlr_despesa_financeira, -vlr_faturamento_bruto,
    -vlr_comissao, -vlr_despesa_venda, -vlr_frete_unitario,
    -nom_cidade, -vlr_faturamento_sem_imposto, -vlr_impostos,
    -vlr_taxa_nota, -vlr_faturamento_sem_ipi, -vlr_frete
  ) %>% 
  filter(cod_uf == UF, cod_site == SITE, cod_tipo_embalagem == EMBALAGEM) %>% 
  select(-cod_uf, -cod_site, -cod_tipo_embalagem) %>% 
  arrange(dat_movimento) %>% 
  pad_by_time(.date_var = dat_movimento, .by = "day", .pad_value = 0) %>% 
  summarise_by_time(.date_var = dat_movimento, .by = "day", 
                    preco_liq_medio = mean(vlr_preco_liquido), 
                    volume_total = sum(qtd_volume)) %>% 
  mutate(
    preco_liq_medio = zoo::na.locf(preco_liq_medio),
    volume_total    = replace_na(volume_total, 0)
  )

skim_without_charts(ibema)

Data summary

Name	ibema
Number of rows	1179
Number of columns	3
_______________________
Column type frequency:
Date	1
numeric	2
________________________
Group variables	None

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
dat_movimento	0	1	2018-01-03	2021-03-26	2019-08-15	1179

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100
preco_liq_medio	0	1	2364.38	1951.63	0	0	3469.03	4027.81	5168.57
volume_total	0	1	23.37	28.98	0	0	13.48	38.13	195.79

Feature Engineering

A etapa anterior resultou em uma série temporal cujos últimos 365 dias são mostrados a seguir, a título de ilustração:

ibema %>% 
  tail(365) %>% 
  plot_time_series(.date_var = dat_movimento, .value = volume_total)

Da análise visual dessa série temporal, algumas características se destacam, a saber:

• Parece haver uma sazonalidade semanal

• Há picos não sazonais cuja explicação pode demandar o cruzamento desses dados com outros dados, por exemplo: eventos de marketing

Para confrontar as hipóteses descritas anteriormente, mostram-se a seguir os gráficos de autocorrelação e de autocorrelação parcial:

ibema %>% 
  plot_acf_diagnostics(.date_var = dat_movimento, 
                       .value = log1p(volume_total))

ibema %>% 
  plot_seasonal_diagnostics(.date_var = dat_movimento, 
                            .value = log1p(volume_total))

Da análise dos gráficos anteriores, nota-se a sazonalidade semanal (7 dias) e também um certo padrão sazonal ao longo dos dias da semana. Essas informações serão importantes no processo de construção dos modelos matemáticos.

Forecasting

calibrate_and_plot <- function(..., type = "testing") {
  
  if (type == "testing") {
    new_data <- testing(splits)
  } else {
    new_data <- training(splits) %>% drop_na()
  }
  
  calibration_tbl <- modeltime_table(...) %>% 
    modeltime_calibrate(new_data)
  
  print(modeltime_accuracy(calibration_tbl))
  
  calibration_tbl %>% 
    modeltime_forecast(
      new_data    = new_data,
      actual_data = data_prepared_tbl
    ) %>% 
    plot_modeltime_forecast(.conf_interval_show = FALSE)
  
}

refit_and_plot <- function(...) {
  
  calibration_tbl <- modeltime_table(...) %>% 
    modeltime_calibrate(testing(splits))
  
  refit_tbl <- calibration_tbl %>%
    modeltime_refit(data_prepared_tbl)
  
  refit_tbl %>%
    modeltime_forecast(new_data    = forecast_tbl,
                       actual_data = data_prepared_tbl) %>%
    
    # Invert Transformation
    mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
      x    = .,
      mean = std_mean,
      sd   = std_sd
    ))) %>%
    mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
      x           = ., 
      limit_lower = limit_lower, 
      limit_upper = limit_upper, 
      offset      = offset
    ))) %>%
    
    plot_modeltime_forecast()
  
}

Transformação dos dados

Em face das análises anteriores, podem ser usados os lags de 7, 14 e 21 dias como tentativas de variáveis adicionais, além do dia da semana e do mês, visando melhorar o desempenho dos modelos matemáticos.

Adicionalmente, a fim de melhorar o desempenho dos modelos matemáticos de machine learning, tomou-se o logaritmo do volume diário, que posteriormente foi normalizado (média zero e desvio padrão igual a um).

# Transformação dos dados
ibema_transformed_tbl <- ibema %>%
  
  # Preprocessamento do Alvo
  mutate(vol_tot_trans = log_interval_vec(volume_total, limit_lower = 0,
                                          offset = 1)) %>%
  mutate(vol_tot_trans = standardize_vec(vol_tot_trans)) %>%
  select(-volume_total, -preco_liq_medio)

# Parâmetros chave
limit_lower <- 0
limit_upper <- 216.3745
offset      <- 1
std_mean    <- -3.17308718389508
std_sd      <- 1.92403817890689

ibema_transformed_tbl %>% glimpse()

## Rows: 1,179
## Columns: 2
## $ dat_movimento <date> 2018-01-03, 2018-01-04, 2018-01-05, 2018-01-06, 2018-01~
## $ vol_tot_trans <dbl> 0.77029164, -1.14306006, -1.14306006, -1.14306006, -1.14~

Visualização do efeito de feature engineering

A fim de testar as hipóteses que motivaram as transformações de dados descritas anteriormente, a seguir mostra-se uma regressão linear do volume (transformado) em função das datas e de:

• lags: 7, 14 e 21 dias

• dia da semana: de segunda-feira a domingo

• mês: de janeiro a dezembro

plot_time_series_regression(
  .data = ibema_transformed_tbl %>% 
    tk_augment_lags(.value = vol_tot_trans, .lags = c(7, 14, 21)) %>% 
    drop_na(),
  .date_var = dat_movimento,
  .formula = vol_tot_trans ~ as.numeric(dat_movimento)
  + wday(dat_movimento, label = TRUE)
  + month(dat_movimento, label = TRUE)
  + .
  - dat_movimento,
  .show_summary = TRUE
)

## 
## Call:
## stats::lm(formula = .formula, data = .data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76844 -0.31520  0.03258  0.48072  2.37820 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            2.866e-01  1.279e+00   0.224 0.822748
## as.numeric(dat_movimento)             -1.593e-05  7.052e-05  -0.226 0.821363
## wday(dat_movimento, label = TRUE).L    2.354e-01  6.210e-02   3.792 0.000158
## wday(dat_movimento, label = TRUE).Q   -1.602e+00  1.004e-01 -15.954  < 2e-16
## wday(dat_movimento, label = TRUE).C    8.591e-02  6.105e-02   1.407 0.159663
## wday(dat_movimento, label = TRUE)^4   -5.819e-01  6.759e-02  -8.609  < 2e-16
## wday(dat_movimento, label = TRUE)^5   -2.290e-03  6.088e-02  -0.038 0.970007
## wday(dat_movimento, label = TRUE)^6   -6.125e-02  6.089e-02  -1.006 0.314676
## month(dat_movimento, label = TRUE).L   7.952e-02  8.026e-02   0.991 0.321998
## month(dat_movimento, label = TRUE).Q   8.877e-02  8.179e-02   1.085 0.277968
## month(dat_movimento, label = TRUE).C  -8.202e-02  8.079e-02  -1.015 0.310188
## month(dat_movimento, label = TRUE)^4  -2.399e-01  8.010e-02  -2.995 0.002802
## month(dat_movimento, label = TRUE)^5   8.740e-02  8.012e-02   1.091 0.275564
## month(dat_movimento, label = TRUE)^6  -4.255e-02  8.018e-02  -0.531 0.595767
## month(dat_movimento, label = TRUE)^7  -1.950e-02  7.912e-02  -0.247 0.805337
## month(dat_movimento, label = TRUE)^8   3.009e-04  7.927e-02   0.004 0.996972
## month(dat_movimento, label = TRUE)^9   2.785e-02  8.097e-02   0.344 0.730918
## month(dat_movimento, label = TRUE)^10 -1.385e-01  8.148e-02  -1.700 0.089418
## month(dat_movimento, label = TRUE)^11  1.390e-01  8.184e-02   1.698 0.089788
## vol_tot_trans_lag7                     2.316e-02  2.986e-02   0.776 0.438083
## vol_tot_trans_lag14                   -3.405e-02  2.989e-02  -1.139 0.254863
## vol_tot_trans_lag21                   -3.311e-02  2.977e-02  -1.112 0.266311
##                                          
## (Intercept)                              
## as.numeric(dat_movimento)                
## wday(dat_movimento, label = TRUE).L   ***
## wday(dat_movimento, label = TRUE).Q   ***
## wday(dat_movimento, label = TRUE).C      
## wday(dat_movimento, label = TRUE)^4   ***
## wday(dat_movimento, label = TRUE)^5      
## wday(dat_movimento, label = TRUE)^6      
## month(dat_movimento, label = TRUE).L     
## month(dat_movimento, label = TRUE).Q     
## month(dat_movimento, label = TRUE).C     
## month(dat_movimento, label = TRUE)^4  ** 
## month(dat_movimento, label = TRUE)^5     
## month(dat_movimento, label = TRUE)^6     
## month(dat_movimento, label = TRUE)^7     
## month(dat_movimento, label = TRUE)^8     
## month(dat_movimento, label = TRUE)^9     
## month(dat_movimento, label = TRUE)^10 .  
## month(dat_movimento, label = TRUE)^11 .  
## vol_tot_trans_lag7                       
## vol_tot_trans_lag14                      
## vol_tot_trans_lag21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7831 on 1136 degrees of freedom
## Multiple R-squared:  0.3994, Adjusted R-squared:  0.3883 
## F-statistic: 35.97 on 21 and 1136 DF,  p-value: < 2.2e-16

Do sumário e do gráfico anteriores, depreende-se que os trabalhos de feature engineering resultaram em um excelente ponto de partida para os algoritmos de machine learning.

Métricas de comparação de desempenho

Daqui em diante, serão testadas nove abordagens para prever os volumes dos próximos 14 dias e os modelos resultantes serão comparados, principalmente, pelas seguintes métricas:

• MAE: Mean Absolute Error

• RMSE: Root Mean Squared Error

• RSQ: R²

Resultados

O gráfico a seguir mostra os dados de treinamento e os de teste (14 dias) usados para treinar e verificar o desempenho dos modelos de forecasting selecionados.

# 1.0 PASSO 1 - CRIAR FULL DATA SET
# - Estender pela janela futura
# - Feature engineering
# - Regressores externos

horizon    <- 14 # 14 dias
lag_period <- 14
rolling_periods <- c(7, 14, 21)
dot_value <- str_glue("vol_tot_trans_lag{lag_period}")

data_prepared_full_tbl <- ibema_transformed_tbl %>%
  
  # Janela futura
  bind_rows(
    future_frame(.data = ., .date_var = dat_movimento, .length_out = horizon)
  ) %>%
  
  # Lags autocorrelacionados
  tk_augment_lags(vol_tot_trans, .lags = lag_period) %>%
  
  # Rolling features to stabilize trend from M5 Competition
  tk_augment_slidify(
    .value   = rlang::sym(dot_value),
    .f       = mean,
    .period  = rolling_periods,
    .align   = "center",
    .partial = TRUE
  ) %>%
  
  # Format Columns
  rename_with(.cols = contains("lag"), .fn = ~ str_c("lag_", .x))

# 2.0 PASSO 2 - SEPARAR os DADOS DE MODELAGEM E DE FORECAST

data_prepared_tbl <- data_prepared_full_tbl %>%
  filter(!is.na(vol_tot_trans))

forecast_tbl <- data_prepared_full_tbl %>%
  filter(is.na(vol_tot_trans))

# 3.0 TRAIN/TEST (DATASET DE MODELAGEM)

splits <- time_series_split(data_prepared_tbl, assess = horizon, 
                            cumulative = TRUE)

splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(dat_movimento, vol_tot_trans)

# 4.0 RECIPES
# - Time Series Signature - Adds bulk time-based features
# - More Transformation to Time Series Signature vars
# - Interactions, if any
# - Fourier Features

recipe_spec_base <- recipe(vol_tot_trans ~ ., data = training(splits)) %>%
  
  # Time Series Signature
  step_timeseries_signature(dat_movimento) %>%
  step_rm(matches("(iso)|(xts)|(hour)|(minute)|(second)|(am.pm)")) %>%
  step_rm(ends_with("(day)|(day7)")) %>% 
  
  # Standardization
  step_normalize(matches("(index.num)|(year)")) %>%
  
  # Dummy Encoding (One Hot Encoding)
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  
  # Fourier
  step_fourier(dat_movimento, period = rolling_periods, K = 2)

# recipe_spec_base %>% prep() %>% juice() %>% glimpse()

# 5.0 LINEAR REGRESSION MODEL

# * Model spec

model_spec_lm <- linear_reg() %>% 
  set_engine("lm")

# * Model Workflow

recipe_spec_lm <- recipe_spec_base %>% 
  step_rm(dat_movimento) %>% 
  step_naomit(contains("lag"))

workflow_fit_lm <- workflow() %>%
  add_model(model_spec_lm) %>%
  add_recipe(recipe_spec_lm) %>%
  fit(training(splits))

# 6.0 GLM REGRESSION MODEL

# * GLM Model Spec

model_spec_glm <- linear_reg(penalty = 0.1, mixture = 0.5) %>%
  set_engine("glmnet")

# * GLM Model Workflow

recipe_spec_glm <- recipe_spec_base %>% 
  step_rm(dat_movimento) %>% 
  step_naomit(contains("lag"))

workflow_fit_glm <- workflow() %>%
  add_model(model_spec_glm) %>%
  add_recipe(recipe_spec_glm) %>%
  fit(training(splits))

# 7.0 RANDOM FOREST MODEL

# * RF Model Spec

model_spec_rf <- rand_forest() %>%
  set_engine("randomForest")

# * RF Model Workflow

workflow_fit_rf <- workflow() %>%
  add_model(model_spec_rf) %>%
  add_recipe(recipe_spec_glm) %>%
  fit(training(splits))

# 8.0 ARIMA MODEL

model_fit_auto_arima <- arima_reg() %>%
  set_engine("auto_arima") %>%
  fit(
    vol_tot_trans ~ dat_movimento 
    + fourier_vec(dat_movimento, period = 7)
    + fourier_vec(dat_movimento, period = 14)
    + fourier_vec(dat_movimento, period = 21),
    data = training(splits)
  )

# 9.0 PROPHET MODEL, NO XREGS

model_fit_prophet <- prophet_reg(
  # changepoint_num    = 25,
  # changepoint_range  = 0.8,
  # seasonality_yearly = FALSE,
  # seasonality_weekly = TRUE,
  # seasonality_daily  = TRUE
) %>%
  set_engine("prophet") %>%
  fit(vol_tot_trans ~ dat_movimento, data = training(splits))

# 10.0 PROPHET MODEL, XREGS

model_fit_prophet_xregs <- prophet_reg(
  # changepoint_num    = 25,
  # changepoint_range  = 0.8,
  # seasonality_yearly = FALSE,
  # seasonality_weekly = TRUE,
  # seasonality_daily  = TRUE
) %>%
  set_engine("prophet") %>%
  fit(
    vol_tot_trans ~ dat_movimento 
    + fourier_vec(dat_movimento, period = 7)
    + fourier_vec(dat_movimento, period = 14)
    + fourier_vec(dat_movimento, period = 21),
    data = training(splits)
  )

# 11.0 CUBIST

model_spec_cubist <- cubist_rules() %>% 
  set_engine("Cubist")

set.seed(123)
workflow_fit_cubist <- workflow() %>% 
  add_model(model_spec_cubist) %>% 
  add_recipe(recipe_spec_lm) %>% 
  fit(training(splits))

# 12.0 XGBOOST

workflow_fit_xgboost <- workflow() %>%
  add_model(
    boost_tree() %>% 
      set_engine("xgboost")
  ) %>%
  add_recipe(recipe_spec_lm) %>%
  fit(training(splits))

# 13.0 NNETAR

model_spec_nnetar <- nnetar_reg(
  # non_seasonal_ar = 1,
  # seasonal_ar     = 3,
  # hidden_units    = 10,
  # penalty         = 10,
  # num_networks    = 20,
  # epochs          = 50
) %>%
  set_engine("nnetar")

set.seed(123)
workflow_fit_nnetar <- workflow() %>%
  add_model(model_spec_nnetar) %>%
  add_recipe(recipe_spec_base) %>%
  fit(training(splits) %>% drop_na())

A tabela a seguir mostra os resultados do desempenho dos modelos matemáticos treinados nos 14 dias de teste, ou seja, nos 14 dias que não participaram do treinamento dos algoritmos.

Os algoritmos comparados a seguir são os seguintes:

• Regressão linear

• GLMNET: modelos lineares generalizados via máxima verossimilhança penalizada

• Random Forest

• ARIMA: média móvel autoregressiva integrada

• PROPHET: algoritmo desenvolvido e utilizado pelo Facebook

• CUBIST: algoritmo baseado em árvores, como o Random Forest, mas com equações lineares nos nós terminais

• XGBOOST: Extreme Gradient Boosting

• NNETAR: uma combinação de redes neurais e de modelos autorregressivos

# 13.0 MODELTIME WORKFLOW

# * Modeltime Table
submodels_tbl <- modeltime_table(
  workflow_fit_lm, workflow_fit_glm, workflow_fit_rf,
  model_fit_auto_arima, model_fit_prophet, 
  model_fit_prophet_xregs, workflow_fit_cubist, 
  workflow_fit_xgboost, workflow_fit_nnetar
)

# * Calibrate Testing Data
submodels_calibrated_tbl <- submodels_tbl %>%
  modeltime_calibrate(testing(splits))

# * Measure Test Accuracy
submodels_calibrated_tbl %>% 
  modeltime_accuracy() %>% 
  table_modeltime_accuracy(
    .interactive = TRUE,
    bordered = TRUE, 
    resizable = TRUE
  )

Modelo Ensemble

Da tabela anterior, excetuando-se o modelo ARIMA, notam-se os excelentes desempenhos dos demais modelos. Uma estratégia para construir um modelo ainda melhor, a partir dos treinados até o momento, consiste em criar um modelo ensemble que calcule suas previsões a partir da média dos três “melhores” modelos, a saber: RANDOMFOREST, XGBOOST e PROPHET.

Os resultados dessa estratégia são mostrados a seguir:

# 15.0 ENSEMBLE

# * Make Ensemble
ensemble_fit_mean <- submodels_tbl %>% 
  filter(.model_desc %in% c("RANDOMFOREST", "XGBOOST", "PROPHET")) %>% 
  ensemble_average(type = "mean")

# * Modeltime Table
ensemble_tbl <- modeltime_table(
  ensemble_fit_mean
)

# * Ensemble Test Accuracy
ensemble_tbl %>%
  combine_modeltime_tables(
    submodels_tbl %>% 
      filter(.model_desc %in% c("RANDOMFOREST", "XGBOOST", "PROPHET"))
  ) %>%
  modeltime_accuracy(testing(splits)) %>% 
  table_modeltime_accuracy(
    .interactive = TRUE,
    bordered = TRUE, 
    resizable = TRUE
  )

Finalmente, os gráficos a seguir mostram, respectivamente, as previsões de volumes nos próximos 14 dias calculadas pelo modelo ensemble, comparadas com os dados de teste, bem como nos próximos 14 dias do horizonte de previsão:

ensemble_tbl %>%
  modeltime_calibrate(testing(splits)) %>% 
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = data_prepared_tbl,
    keep_data   = TRUE
  ) %>%
  # Invert Transformation
  mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
    x    = .,
    mean = std_mean,
    sd   = std_sd
  ))) %>%
  mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
    x           = ., 
    limit_lower = limit_lower, 
    limit_upper = limit_upper, 
    offset      = offset
  ))) %>%
  plot_modeltime_forecast()

ensemble_refit_tbl <- ensemble_tbl %>%
  modeltime_calibrate(testing(splits)) %>% 
  modeltime_refit(data_prepared_tbl)

# * Visualize Ensemble Forecast
ensemble_refit_tbl %>%
  
  modeltime_forecast(
    new_data    = forecast_tbl,
    actual_data = data_prepared_tbl,
    keep_data   = TRUE
  ) %>%
  
  # Invert Transformation
  mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
    x    = .,
    mean = std_mean,
    sd   = std_sd
  ))) %>%
  mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
    x           = ., 
    limit_lower = limit_lower, 
    limit_upper = limit_upper, 
    offset      = offset
  ))) -> fore_vol
  
plot_modeltime_forecast(fore_vol)

Forecasting dos Preços Líquidos

Como o objetivo desta PoC é prever o faturamento líquido nos próximos 14 dias - em um estado, site e para uma determinada embalagem - , nos faltam as previsões correspondentes dos preços líquidos.

A bem da concisão, não repetiremos todo o processo anterior para os preços líquidos, cujos dados passaram por processo semelhante ao apresentado neste relatório para os volumes diários.

Modelou-se a evolução dos preços líquidos com as mesmas técnicas e algoritmos utilizados para os volumes, cujos resultados constam na tabela a seguir:

precos <- read_rds("00_artifacts/price_artifacts.rds")

# * Calibrate Testing Data ----
submodels_calibrated_tbl <- precos$submodels_tbl %>%
  modeltime_calibrate(precos$test_tbl)

# * Measure Test Accuracy ----
submodels_calibrated_tbl %>% 
  modeltime_accuracy() %>% 
  table_modeltime_accuracy(
    .interactive = TRUE,
    bordered = TRUE, 
    resizable = TRUE
  )

Da análise da tabela anterior, depreende-se que o melhor modelo é o NNETAR, sendo que os gráficos a seguir mostram, respectivamente, as previsões de preços líquidos nos próximos 14 dias calculadas pelo modelo selecionado, comparadas com os dados de teste, bem como nos próximos 14 dias do horizonte de previsão:

best_model_fit <- precos$best_model_fit

calibration_tbl <- modeltime_table(best_model_fit) %>% 
  modeltime_calibrate(new_data = precos$test_tbl)

calibration_tbl %>% 
  modeltime_forecast(
    new_data    = precos$test_tbl,
    actual_data = precos$data_prepared_tbl
  ) %>% 
  # Invert Transformation
  mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
    x    = .,
    mean = precos$inv_params$std_mean,
    sd   = precos$inv_params$std_sd
  ))) %>%
  mutate(across(.value:.conf_hi, .fns = exp
  )) %>%
  
  plot_modeltime_forecast()

refit_tbl <- calibration_tbl %>%
  modeltime_refit(precos$data_prepared_tbl)

# Visualize Forecast
refit_tbl %>%
  
  modeltime_forecast(
    new_data    = precos$fore_tbl,
    actual_data = precos$data_prepared_tbl,
    keep_data   = TRUE
  ) %>%
  
  # Invert Transformation
  mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
    x    = .,
    mean = precos$inv_params$std_mean,
    sd   = precos$inv_params$std_sd
  ))) %>%
  mutate(across(.value:.conf_hi, .fns = exp
  )) -> fore_price
  
plot_modeltime_forecast(fore_price)

Forecasting do Faturamento Líquido

De posse das previsões de volumes e de preços, calculam-se as previsões de faturamento por meio da multiplicação de ambas, conforme mostrado no gráfico a seguir:

fore_vol %>% 
  tail(180) %>% 
  # dplyr::filter(.key == "prediction") %>% 
  select(dat_movimento, .key, .value) %>% 
  rename(vol = .value) %>% 
  bind_cols(
    fore_price %>% 
      tail(180) %>% 
      # dplyr::filter(.key == "prediction") %>% 
      select(.value) %>% 
      rename(price = .value)
  ) %>% 
  mutate(net_income = vol * price) %>% 
  pivot_longer(cols = c(vol, price, net_income)) %>% 
  plot_time_series(.date_var = dat_movimento, .value = value, 
                   .color_var = .key, .facet_vars = name, .smooth = FALSE)

Conclusões e Recomendações

A aplicação das técnicas mencionadas anteriormente para os dados do estado de São Paulo, bobinas do site TVO, permite concluir que os modelos preditivos desenvolvidos apresentaram muito bom desempenho e que, portanto, podem ser usados para prever acuradamente tanto o preço líquido quanto o volume de vendas em um horizonte de até 14 dias. Essas previsões, se estendidas para outras regiões, sites e tipos de embalagem, certamente poderão apoiar o processo decisório da Ibema no que diz respeito ao planejamento e controle da produção, bem como em relação às estratégias de marketing.

Em face das conclusões anteriores, recomenda-se o seguinte:

• Avaliar modelos para outras unidades da federação, tipos de embalagem, cidades, etc.

• Listar os forecastings por ordem decrescente de faturamento líquido previsto, conforme a granularidade que se deseje avaliar, e.g. estado, cidade, tipo de embalagem, configuração de produto, etc.

• Automatizar o processo por meio de um sistema HPFS - High Performance Forecasting System